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We present results for light meson masses and pseudoscalar decay constants from the first of a 
series of lattice calculations with 2+1 dynamical flavors of domain wall fermions and the Iwasaki 
gauge action. The work reported here was done at a fixed lattice spacing of about 0.12 fm on 
a 16 3 x 32 lattice, which amounts to a spatial volume of (2 fm) 3 in physical units. The number 
of sites in the fifth dimension is 16, which gives m rC s = 0.00308(4) in these simulations. Three 
values of input light sea quark masses, m; ca m 0.85 m s , 0.59 m s and 0.33 m 3 were used to allow 
for extrapolations to the physical light quark limit, whilst the heavier sea quark mass was fixed 
to approximately the physical strange quark mass m s . The exact rational hybrid Monte Carlo 
algorithm was used to evaluate the fractional powers of the fermion determinants in the ensemble 
generation. We have found that f„ = 127(4) MeV, f K = 157(5) MeV and f K /f„ = 1.24(2), where 
the errors are statistical only, which are in good agreement with the experimental values. 

PACS numbers: 11.15.Ha, 11.30.Rd, 12. 38. Aw, 12.38.-t 12.38.Gc 



I. INTRODUCTION 

The successful determination of the light hadron spec- 
trum serves as an important test of lattice QCD. With 
the emergence of powerful computers and new algo- 
rithms, lattice QCD has entered an era of dynamical sim- 
ulations where two light quarks and one strange quark 
are included in the vacuum polarization effects. (These 
simulations are referred to as 2+1 flavor simulations.) 
Since current simulations use light quarks that are heav- 
ier than the physical values, extrapolations to the light 
quark region are still required, and having control of these 
extrapolations is vital for a comparison with experimen- 
tal results. Chiral perturbation theory provides a frame- 
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work for such extrapolations, although the range of quark 
masses for which a given order of chiral perturbation the- 
ory is accurate is still under investigation. 

Symmetry plays an important role in hadron physics. 
The domain wall fermion (DWF) formulation [H, 0,3 re- 
spects flavor symmetry and has approximate chiral sym- 
metry, with the introduction of an auxiliary fifth dimen- 
sion. When the extent of this fifth dimension, denoted 
as L s , goes to infinity, chiral symmetry is fully recov- 
ered. In addition, unlike Wilson and staggered fermions, 
domain wall fermions (and the closely connected overlap 
fermions) allow chiral symmetry to be recovered at finite 
lattice spacing. While certainly an important improve- 
ment for the light hadron spectrum calculation presented 
here, this good chiral symmetry is vital for accurate oper- 
ator renormalization and control of operator mixing in a 
wide variety of important hadronic matrix elements. As 
evidenced by existing work in quenched QCD, the chi- 
ral properties of domain wall fermions make calculations 
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possible that cannot currently be done with other formu- 
lations 0, [j| . To move beyond the quenched approxima- 
tion, 2+1 flavor ensembles are needed. The production 
of these ensembles and the determination of their basic 
properties are important steps towards the goal of mea- 
suring a wide variety of physically interesting quantities 
using this formulation and are the focus of this paper. 

In practice, L s is taken to be on the order of 10-20 
sites so that some residual chiral symmetry breaking re- 
mains. However, this residual chiral symmetry break- 
ing is highly suppressed. Compared to Wilson fermions 
where the residual chiral symmetry breaking is an 0(1) 
effect, for DWF it is of C(10 -3 ) or smaller, depending 
on the size of L s . To leading order in the lattice spac- 
ing a, the only effect of this residual chiral symmetry 
breaking is a simple mass renormalization (denoted by 
TO rC s) added to the input quark mass. Additional chiral 
symmetry breaking effects enter at higher order in a, but 
they are suppressed by the additional powers of a and 
also vanish as L s becomes large. This situation simpli- 
fies the chiral extrapolation of physical observables, and 
allows us to study those aspects of hadron phenomenol- 
ogy where chiral symmetry is important with reduced 
systematic uncertainties. 

Domain wall fermions arc both on- and off-shell im- 
proved, since for any quantity a traversal of the fifth 
dimension is required to mix chiralities. Such a traver- 
sal is suppressed for large L s , a suppression which gives 
rise to the small value for to 1cs and which depends on 
L s as exp(— aL s ) if the generic, weak coupling behav- 
ior of residual chiral symmetry breaking effects is used 
Q. Thus, the resulting theory is accurate up to terms of 
0(aexp(-(xL s )) and 0(a 2 ). 

For amplitudes which violate chiral symmetry by two 
units, a suppression of order m^ es or exp(— 2aL s ) is ex- 
pected. This estimate may fail in the case of compos- 
ite operators with many fields at the same space-time 
point fQj, where one must give special consideration to chi- 
ral symmetry breaking arising from rare, localized modes 
that are undamped in the fifth dimension. However, the 
detailed flavor structure of the operator must be consid- 
ered which, in important cases, prevents these localized 
modes from contributing Q. The result is that chirally 
disallowed operator mixings are generally under good 
control with DWF at finite L s . 

For power divergent quantities, one finds that if an op- 
erator receives an mf/a n contribution, where to,/ is the 
quark mass, then for finite L s it will generally receive 
an 0(m Ies )/a n contribution Because the residual 
chiral symmetry breaking is entering a divergent ampli- 
tude, it will not be determined by the chiral symmetry 
breaking at low energies described by TO res - However, 
it will be of this order, as has been demonstrated in 
[||. Thus, for such quantities the simple replacement 
mf/a n — > (to/ + m rcs )/a™ provides a sensible estimate of 
but not an accurate value for the effects of residual sym- 
metry breaking on the quantity in question. For exam- 
ple, such a term enters the chiral condensate (qq). This 



implies that to compute (qq) with domain wall fermions 
an extrapolation to large L s must be made or that this 
physical quantity should be determined from the den- 
sity of Dirac eigenvalues at zero, so that ultraviolet chi- 
ral symmetry breaking effects do not enter. In contrast, 
a similar, uncertain 0(TO res /a 2 ) term which appears in 
the matrix elements of the weak interaction operator Og , 
does not effect the quantity of interest which is propor- 
tional to to/ For a more complete discussion of the 
effects of residual chiral symmetry breaking for domain 
wall fermion see, for example, Refs 0, El S> 1| • 

These benefits of controlled chiral symmetry breaking 
must be weighed against the computing cost for domain 
wall fermions, which are naively 0(L S ) more expensive 
than conventional fermion formulations [101 ]. In balance, 
we believe that the benefits of domain wall fermions out- 
weigh their extra numerical costs, because of the much 
larger range of observables that are accessible to this for- 
mulation. 

Studies which established a set of parameter values 
suitable for 2+1 flavor simulations with domain wall 
fermions were reported in [111 ]. In the present paper we 
describe the first of a series of 2+1 flavor DWF simu- 
lations which use the parameters as determined in [ill ], 
and explore in detail the systematic effects arising in such 
a full dynamical simulation with domain wall fermions. 
Specifically, we use the Iwasaki gauge action with gauge 
coupling (3 — 2.13 on a lattice of 16 3 x 32, which gives a 
lattice spacing of about 0.12 fm and a spatial lattice ex- 
tent of about 2 fm. This volume should be large enough 
to give accurate results for the physics of the pseudoscalar 
mesons simulated here: those with a mass ratio to vec- 
tor meson masses in the range 0.5 < mp/my < 1.0. 
Larger volume simulations at the same coupling using 
lighter quark masses are underway, as are simulations at 
weaker couplings, also with lighter quark masses. The 
volumes used in the work presented here do not allow 
lighter dynamical quark masses to be used, so the full 
benefits of DWF will not be visible. However, the cur- 
rent simulations allow us to demonstrate that more costly 
simulations at larger volumes and lighter quark masses 
are feasible. Since baryonic observables may suffer from 
non-negligible finite volume effects in the current volume, 
we will only focus on light meson physics in this paper. In 
particular, we will present results for light pseudoscalar 
meson masses and decay constants, and calculate quanti- 
ties of direct phenomenological interest, such as the light 
quark masses and fx/ fir- 

One of the many complexities of 2+1 flavor simulations 
is the presence of the fractional powers of the fermion de- 
terminants in the path integral, which the conventional 
hybrid Monte Carlo algorithm is not able to evaluate. We 
employed the rational hybrid Monte Carlo (RHMC) al- 
gorithm by Clark and Kennedy [H, [H, 03] , which is not 
only free of finite step size errors, but also has proved 
to be quite efficient after recent improvements and tun- 
ings [15j. The gauge configurations used in this work 
were mainly generated with one variant of the RHMC 
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algorithm. We describe the details of this algorithm in 
Section |IIJ The algorithm was further improved after the 
major part of this work had been done. It is thus to our 
advantage to study how the computing cost and ensem- 
ble properties change with the new improvements. We 
compare these two variants of the RHMC algorithm in 
Section ED 

This paper is organized as follows. In Section [IT] we 
present the details of generation of the gauge configura- 
tions using the RHMC algorithm. We discuss the ther- 
malization and autocorrelations of the ensembles in Sec- 
tion [TTTJ Section HVl gives the results of the light meson 
masses and decay constants, including the residual chiral 
symmetry breaking in these simulations. The results for 
physical observables in the chiral limit are presented in 
Section [Vj Section I VII discuses an extension of one of 
the ensembles using an improved RHMC algorithm. Our 
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conclusions are given in Section IVlIl 



II. GENERATION OF GAUGE 
CONFIGURATIONS 



In this study we have used the Iwasaki 
renormalization-group improved gauge action and 
2 + 1 flavors of dynamical domain wall fermions; details 
of the action and our notation may be found in [ill ]. 
Our ensembles were generated using two new variants of 
the RHMC algorithm as described below. 

For the case of two degenerate light quarks of mass, 
to/, and a strange quark of mass, m s , integrating out the 
fermion fields in the path integral leads to 



det 


D dwf ( M 5,mi) D D wf {M 5 , TO/ ) 


det 1 / 2 


d dwf (M 5 ,m s ) Ddwf (M 5 ,m s ) 


det 


4wf( m 5, 1)£W(M 8 , 1) 


det 1 / 2 


D^ WF (M 5 ,l)D DWF (M 5 ,l) 
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where -Ddwf (-^5, m) is the domain wall Dirac opera- 
tor, M5 is the height of the domain wall, and m is the 
fermion mass. The determinants in the denominator 
arise from the Pauli-Villars fields introduced to cancel 
the bulk infinity which arises as the size of the fifth di- 
mension L s — > 00. 

We may use the identity det(M) = [det 1/n (M)] n 
to simulate various different decompositions of 
the same fermionic determinants using RHMC. 
If we adopt the shorthand notation D(mj) = 
Dl )WF (M 5 ,m i )D I ) WF (M 5 ,m i ), we may write the 
ratio of determinants in Eq. ([1]) as 



det 1/2 [P(TO/)] det 1/2 [P(TO/)] det 1/2 [P(TOs)] 
det[P(l)] det 1/2 [P(l)] 



(2) 



In our previous paper [111 ] , each determinant in Eq. @ 
was associated with a separate pseudo-fermion field, and 
a standard RHMC algorithm 0, [l7[ was employed. 



The pseudo-fermion action for molecular dynamics 
evolution of a single quark flavor of mass to takes the 
form S p { — 0^P _1 / 2 (m)(/). We may improve the effi- 
ciency of RHMC by taking advantage of the properties 
of the expansion coefficients in the partial-fraction formu- 
lation of the p th -order rational approximation, rr^ 1 / 2 ~ 

|ao + J2k=i x+f3 k }• T ne "shifts", fa, which are all real 



and positive, are ordered according to increasing magni- 
tude, and the ak are all real and positive. The fermion 
force is 



S' pf =J2 a k^ (P(m) + faY 1 V'{m) (V(m) + fay 1 <j> . 

I 



(3) 



fc=i 



where the prime denotes differentiation with respect to 
the gauge field. Somewhat remarkably, the inverses with 
the smallest shifts, fa, which are the most expensive 
to evaluate, are associated with the smallest residues 
ak, and thus lead to relatively-small contributions to 
the fermion force [3, E3]- ft 1S then natural to split 



the sum in Eq. ([3]) into two parts: "light poles" with 
1 < k < r, and "heavy poles" with r < k < p, and eval- 
uate their contributions to the fermion force on different 
timescales using a Sexton- Weingarten multi-timescale in- 
tegrator [l8j]. The light-pole contributions are evaluated 
only on the largest timescale, the heavy poles on an in- 
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termediate timescale, and the gauge fields on the finest 
timescale. 



TABLE I: Parameter values for 2 + 1 flavor RHMC runs on 
16 3 x 32 lattices with L s = 16, Ms = 1.8, /3 = 2.13, mf a = 
0.04. Nmu units is the number of molecular dynamics time 
units obtained from an ordered start, except for the mf* = 
0.03 ensemble generated using RHMC II, which started from 
the last configuration generated using RHMC I. N atep i is the 
number of coarse integration steps of size St per MD time 
unit. The intermediate integration scale is 5r/7Vstep2 and the 
finest scale is St / (N ate p2 N a t ep 3). The RHMC acceptance rate 
is given in the second last column. The final column denotes 
the variant of the RHMC algorithm used. 
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RHMC I 


0.03 
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0.14 


7 


2 


8 
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RHMC I 
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3580 


0.25 


4 


1 


6 


80% 


RHMC II 



In order to maximize the acceptance rate per unit 
computational cost, we optimized the order of the ra- 
tional approximations and the points, r, at which the 
partial-fraction sums are split, whilst also varying the 
three Sexton- Weingarten timescales. The optimal ratio- 
nal approximations and splits were (p, r) = (10, 4) for the 



light quarks, (p, r) = (10,3) for the strange quark, and 
(p, r) = (6,1) for the Pauli-Villars field associated with 
the strange quark. 

In the RHMC accept/reject step, we used rational 
approximations of order 16 and 9, respectively, for the 
quark and Pauli-Villars fields. 

For all ensembles, we used the second-order minimum- 
norm (2MN) integrator [H, G]J H3, EH with approxi- 
mately unit-length trajectories. We tuned the free pa- 
rameter in 2MN to A = 0.22 to minimize the RMS change 
in the Hamiltonian, {(AH) 2 } 1 / 2 , along a set of RHMC 
trajectories. 

This completes the description of the RHMC I algo- 
rithm used to generate the ensembles in the main part 
of this work. We used three values of to/ to generate 
the gauge configurations using this algorithm, with m s 
fixed to 0.04 (in lattice units). Further parameters and 
simulation details are given in Table |TJ 

Recently, we developed a more efficient implementa- 
tion, RHMC II, which was used to generate the additional 
ensemble discussed in Section I VII Following Aoki et al 
[22J one may employ a single pseudo-fermion field to esti- 
mate directly each ratio of determinants in Eq. (fTJ) . This 
reduces the stochastic noise in the molecular dynamics 
evolution and speeds up the calculation because a larger 
step size can be used in the integration. Furthermore, we 
can rewrite the ratio of fcrmion determinants as 
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= det 



T>(mi) 
V{m s ) 



det 



2?(1) 



1/2 



det 



V(m s ) 
V{\) 



1/2 



det 



P(m.) 
V(l) 



1/2 



(4) 



thus preconditioning the light quark Dirac operator by 
the strange quark [2l|, [23| . We associate the first determi- 
nant on the RHS of Eq. Q with a single pseudo-fermion 
field, and introduce separate pseudo-fermion fields for 
each of the three "strange quarks" in the rest of the deter- 
minants on the RHS. Again we used a multiple-timescale 
integrator, with the preconditioned light quarks using the 
largest step size, the strange/Pauli-Villars fields an in- 
termediate size, and the gauge fields the finest size. The 
optimal free parameter of the MN2 integration scheme 
was found to be A = 0.22 in this case. The relevant 
parameters are also given in Tabled] 

To distinguish the valence quark masses from the sea 
quark masses, we will use TO| ca to denote the light dynam- 
ical quark masses, and toi and mi to denote the masses 
of two valence quarks that make up a meson. Whenever 
a meson is composed of two degenerate quarks, a short- 
hand notation m va i may be used to refer to the mass 
of each component quark. We use the lattice units for 
various quantities throughout the paper, unless physical 
units are explicitly given. 



III. THERM ALIZATION AND 
AUTOCORRELATIONS 



In this section we will discuss the thermalization and 
autocorrelations of the ensembles which build a founda- 
tion for the data analysis of various quantities. For most 
of our measurements, we discarded the first 500 molecu- 
lar dynamics time units to account for the thermalization 
of the gauge fields. Shown in Figure[T]is the evolution his- 
tory of the plaquette for all three ensembles. The straight 
lines are the average values of the plaquette for each en- 
semble from time unit 2000 to 4000. It is evident that 
the gauge fields have come to equilibrium after O(500) 
molecular dynamics time units. We have also measured 
the chiral condensate (qq) through a stochastic estimator 
with one single hit for the first 3000 molecular dynamics 
time units, starting from time unit 500, and found no 
signs of further equilibration (Figure [2]) . This supports 
our choice of cut to allow for the thermalization of the 
ensemble. 

To obtain reliable estimates for the statistical errors 
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FIG. 1: Evolution history of the plaquette for all three ensem- 
bles, measured every MD time unit. The straight lines are the 
averages from trajectory 2000 to 4000 for each ensemble. 



in the physical observables, we also investigate the auto- 
correlations in the ensembles. While the autocorrelation 
time may differ from quantity to quantity, we calculate 
the autocorrelations for the quantities of our direct inter- 
est, such as meson correlation functions. We adopted the 
same method as in |22| and calculated the integrated au- 
tocorrelation time, ri nt , for the two-point pseudoscalar 
correlation function at time slice 12 obtained from a 
Coulomb gauge fixed wall source at t = with two degen- 
erate valence quarks with m va \ = 0.01 for the mf a = 0.01 
ensemble, where the correlation function was measured 
every 5 molecular dynamics time units. This is shown 
in Figure [3] We can see that the integrated autocor- 
relation time reaches a plateau of about 20 to 25 time 
units. The separation between two independent mea- 
surements is then about 2rj n t ~ 50 molecular dynamics 
time units. Thus we chose to average the data over a 
block of measurements in such a way that the span of 
the measurements in each block covers at least 50 molec- 
ular dynamics time units whenever possible prior to the 
standard jackknife analysis. 

The evolution of the topological charge is another in- 
dicator of how fast the gauge fields decorrelate. On each 
ensemble we measured the topological charge using the 
classically Q(a 2 ) improved definition of the field tensor as 
defined in Figure 2] shows that the gauge field moves 
frequently between different topological charge sectors. 
Also shown in Figure [4] are the histograms of the topo- 
logical charge, which indicate that different topological 
sectors are sampled reasonably well. However, it is also 
worth noting that correlations on a scale of a few hun- 



FIG. 2: Evolution history of the quark condensate (qq) for the 
first 3000 molecular dynamics time units at three dynamical 
points: m va i = m? 63 = 0.01 (top), m va i = = 0.02 (mid- 
dle) and m va i = mf 11 = 0.03 (bottom). The measurements 
were taken every 5 molecular dynamics time units. 
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FIG. 3: Integrated auto-correlation time for the pseudoscalar 
two-point correlation function at time slice 12 with mf a = 
m va i = 0.01. 



dred molecular dynamics time units are present in the 
TO sca _ o Q3 ensemble. This may imply that the value 
for the autocorrelation time discussed in the previous 
paragraph is under-estimated. Ideally we should use a 
larger block size prior to the jackknife analysis to have 
more robust error estimates, but we are constrained by 
the moderate lengths of the simulations and have chosen 
to use a block size that is not too small, but at the same 
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FIG. 4: The topological charge as measured by gluonic fields 
on the three ensembles. The left three panels are the topologi- 
cal charge as a function of the molecular dynamics time. Each 
of the three graphs on the right represents the normalized his- 
togram of the corresponding topological charge evolution on 
the left. 



time leaves a reasonable number of jackknife blocks to 
perform covariant fits when necessary. 



IV. HADRONIC OBSERVABLES 
A. Calculation Details and Fitting Procedures 

We used several different source and sink operators to 
calculate the quark propagators needed for the two-point 
correlation functions. To be specific, we used a local 
operator (denoted as L), a Coulomb gauge fixed wall op- 
erator (denoted as W), and a Gaussian smeared opera- 
tor [U (denoted as G), for the source and/or sink, in an 
attempt to investigate the systematic uncertainties com- 
ing from excited-state contamination. Here we adopt the 
same notation as in [TTI | and denote the types of me- 
son correlation functions by the different source and sink 
operators used in the quark propagators. For example, 
GL-GL means a correlation function composed of two 
quark propagators, each of which is computed using the 
Gaussian smeared source operator and the local sink op- 
erator. The sources are placed at multiple locations to 
reduce the short-term noise induced by the fluctuations 
in the gauge fields. The correlation functions are aver- 
aged over different source locations before the standard 
jackknife analysis is performed. Details of the calculation 
parameters are tabulated in Tables |TT] and IIII1 

In parallel to the spectrum calculations, we have also 
done another set of measurements necessary for the 
extraction of weak matrix elements (see. for example, 
Ref. 3), in which the quark propagators were calculated 
from Coulomb gauge fixed wall sources at time slices 



isrc = 5 and t src = 27. In contrast to the standard spec- 
trum measurement, here the meson correlation functions 
were constructed from the sum of propagators computed 
with periodic (P) and anti-periodic (A) boundary condi- 
tions in the time direction. This has the effect of doubling 
the temporal extent of the lattice which gives a longer 
time range to extract the meson masses. Five different 
valence masses are used to compute the quark propaga- 
tors and the meson correlation functions are constructed 
from all possible combinations of these propagators. The 
details of measurements are given in Table IIVI We will 
use "P+A" to refer to the results obtained from these 
data sets hereafter. 

To get the meson masses and the corresponding am- 
plitudes of the correlation functions, we performed the 
standard covariant x 2 fits to the correlation functions 
using the hyperbolic cosine form: 



C(t) = A 



-m G t , e — ma{T—t) 



(5) 



where ma is the ground-state meson mass, and t is the 
time slice relative to the source. Note that here T is 
the extent in the time direction of the lattice, i.e. 32 
for the standard spectrum measurements and 64 for the 
P+A measurements. This relation assumes that the sep- 
aration between the source and sink operators is large 
enough that the contamination from the excited states is 
negligible. However, different source operators have dif- 
ferent degrees of overlap with excited states, and hence 
the minimum time slice that should be included in each fit 
may differ. To reduce the uncertainties in fit parameters 
from the choice of fit range, we also performed simulta- 
neous fits to two types of correlation functions using a 
double-cosh form which give both the ground-state and 
the first-excited-state meson masses. We found that the 
ground-state meson masses obtained from the double- 



TABLE II: Measurement parameters for the WL-WL corre- 
lation functions for all the three ensembles. mf a denotes the 
light sea quark mass, mi and m,2 denote the two valence 
quarks composing the mesons. A is the separation between 
measurements in molecular dynamics time units. 
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5 


710 


0, 16 


0.03 


0.01 


0.01 


500 - 


3995 


5 


700 


0, 16 


0.03 


0.02 


0.02 


500 - 


3995 


5 


700 


0, 16 


0.03 


0.03 


0.03 


500 - 


3995 


5 


700 


0, 16 


0.03 


0.04 


0.04 


500 - 


3995 


5 


700 


0, 16 
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TABLE III: Measurement parameters for the LL-LL, GG-GG 
and GL-GL correlation functions for all the three ensembles. 
m.; ea denotes the light sea quark mass, mi and m,2 denote the 
two valence quarks composing the mesons. A is the separation 
between measurements in molecular dynamics time units. 





mi 


m.2 


range 


A 


■^V mcas 


t 


src 


0.01 


0.01 


0.01 


500 - 


3990 


10 


350 


0, 


16 


0.01 


0.01 


0.04 


500 - 


3990 


10 


350 


o. 


16 


0.01 


0.04 


0.04 


500 - 


3990 


10 


350 


o. 


16 


0.01 


0.01 


0.01 


505 - 


3995 


10 


350 


8, 


24 


0.01 


0.01 


0.04 


505 - 


3995 


10 


350 


8, 


24 


0.01 


0.04 


0.04 


505 - 


3995 


10 


350 


8, 


24 


0.02 


0.02 


0.02 


500 - 


3990 


10 


350 


0, 


16 


0.02 


0.02 


0.04 


500 - 


3990 


10 


350 


0. 


16 


0.02 


0.04 


0.04 


500 - 


3990 


10 


350 


o. 


16 


0.02 


0.02 


0.02 


505 - 


3995 


10 


350 


8, 


24 


0.02 


0.02 


0.04 


505 - 


3995 


10 


350 


8, 


24 


0.02 


0.04 


0.04 


505 - 


3995 


10 


350 


8, 


24 


0.03 


0.03 


0.03 


500 - 


3990 


10 


350 


0, 


16 


0.03 


0.03 


0.04 


500 - 


3990 


10 


350 


o, 


16 


0.03 


0.04 


0.04 


500 - 


3990 


10 


350 


o, 


16 



TABLE IV: Measurement parameters for the correlation func- 
tions measured with P+A boundary conditions for all the 
three ensembles. m,; ca denotes the light sea quark mass, mi 
and rri2 denote the two valence quarks composing the mesons. 
A is the separation between measurements in molecular dy- 
namics time units. 



mr 


mi 


m,2 


range 


A 


N; mcas 


^src 


0.01 


0.01 


- 0.05 


0.01 


- 0.05 


1000 - 4000 


20 


150 


5. 


27 


0.02 


0.01 


- 0.05 


0.01 


- 0.05 


1000 - 4000 


20 


150 


5. 


27 


0.03 


0.01 


- 0.05 


0.01 


- 0.05 


1000 - 4000 


20 


150 


5. 


27 



cosh fits were consistent with the simple cosh fit in Eq. ([5]) . 
Since we are mostly interested in the ground-state meson 
masses in this paper, we will only present results obtained 
from the single cosh fits unless otherwise stated. 

In order to determine a proper fit range for the meson 
masses using Eq. ([SJ , we examined the effective masses as 
defined by 



m c ft = cosh 



C(t + 1) + C(t-1) 
2C(t) 



(0) 



and chose the minimum time slice in the fit to be the 
onset of the plateau in the effective mass. A more strin- 
gent way to determine the best fit range is to check how 
the mass obtained from the fit and the resulting x 2 /dof 
vary with the fit range. Since the fits are insensitive to 
the maximum time slice used, we only investigated the 
variations with respect to the lower bound (denoted as 
t min ) and fixed the upper bound (denoted as t max ) to T/2 
for the GL-GL and WL-WL data sets. In Figure the 
fitted masses and % 2 /dof for the pseudoscalar correlation 



functions are plotted for each light dynamical point with 
m va i = mf ea . As we can see, the central values of the 
fitted masses are insensitive to the choice of i m i n when 
£ m in is above 8 or 9 for both GL-GL and WL-WL data 
sets. To obtain high confidence level for the fits, we also 
seek to have minimal x 2 /dof in the time range where the 
fits are stable. 



m, = 0.01 



0.252 - 



-P 1 0.248 



0.244 - 



0.24 n L 



o 4 - 







0.33 



0.32- 



■g 6; 

•° 4- 



E "0.39- 



() O 
[] □ 



o Q 



O GL-GL 
□ WL-WL 



10 



10 



12 



[] 

G 



12 



14 



14 



mi sea = 0.02 




ii 



A GL-GL, t src = 8, 24 
+ GL-GL, t = 0, 16 

sre 



O GL-GL, tj c = (),' 8, 16,24 
□ WL-WL, t = 0, 16 



14 



10 



12 



14 



sea 

ir^ = 0.03 




FIG. 5: Variation of the pseudoscalar masses obtained from 
the correlated fits and the corresponding x 2 /dof with respect 
to tmin for each of the dynamical points. 
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Typically the fitted pseudo-scalar meson masses from 
WL-WL and GL-GL correlation functions agree quite 
well when t mm is above 9. An exception is the mp a = 
0.02 data set, where results from GL-GL correlation func- 
tions are larger than those from WL-WL correlation func- 
tions by a few standard deviations. The deviations result 
from the fact that the GL-GL measurements were per- 
formed at an additional set of source locations which were 
not used in the WL-WL measurements. When only the 
data with the source operator at time slices and 16 are 
included in the analysis, the results obtained from GL- 
GL and WL-WL agree within statistical uncertainties, 
whilst the GL-GL measurements from sources at time 
slices 8 and 24 give much larger values. Hence when all 
the sources are combined in the analysis, the values of the 
fitted masses are much higher from the GL-GL correla- 
tors than those from the WL-WL correlators. The fact 
that the results are so sensitive to the different source lo- 
cations indicates that we have not sampled phase space 
extremely well. Thus, the errors obtained from the stan- 
dard jackknife procedure may be under-estimated. 



B. Meson Masses 

The pseudoscalar meson masses can be obtained from 
several different correlation functions [1, We will 

show the results from both the pseudoscalar (PP) and 
axial vector (AA) channels. The effective masses in the 
pseudoscalar channel at the light dynamical points for 
all three ensembles are shown in Figure [Sj The effective 
masses from LL-LL correlation functions typically have 
later approaches to plateaus and leave us with fewer time 
slices to perform the fits. Thus we will not quote results 
from these data sets unless otherwise noted. We also 
note that there are some variations in the effective masses 
greater than estimates of the standard deviation from the 
variance of the data. The presence of these variations 
with 4000 MD time units may be another indication of 
long-term autocorrelations in these ensembles. Although 
we do not fully understand the cause of this problem, it 
is likely due to the algorithm we used to generate these 
gauge configurations. This will be further addressed in 
Section ED 

We determined the vector meson masses by averaging 
over the three polarizations to reduce statistical fluctua- 
tions. Some representative effective masses for the vector 
mesons are shown in Figure [7] They are typically nois- 
ier than the pseudoscalar effective masses, and approach 
plateaus later. As the quark masses become lighter, the 
noise on the effective masses increases. 

As discussed already, we have chosen to fit the corre- 
lation functions with different source/sink operators in- 
dependently to obtain the corresponding meson masses, 
which are given in Appendix [X] The sources differ not 
only in their spatial characteristics, but also in their tem- 
poral location in the lattices. As such, they sample dif- 
ferent parts of the gauge fields. The different combina- 



TABLE V: Results of the pseudoscalar and vector meson 
masses averaged from different measurements. The quoted 
errors are the largest statistical errors of the different mea- 
surements multiplied by 1.5 to account for possible underes- 
timation of errors as described in the text. 



mr 


mi 


7T12 


mp 


mv 


0.01 


0.01 


0.01 


0.247(3) 


0.549(20) 


0.01 


0.02 


0.01 


0.290(3) 


0.564(14) 


0.01 


0.02 


0.02 


0.323(3) 


0.577(11) 


0.01 


0.03 


0.01 


0.325(3) 


0.580(11) 


0.01 


0.03 


0.02 


0.357(3) 


0.595(9) 


0.01 


0.03 


0.03 


0.385(3) 


0.609(8) 


0.01 


0.04 


0.01 


0.356(3) 


0.599(10) 


0.01 


0.04 


0.02 


0.387(3) 


0.611(8) 


0.01 


0.04 


0.03 


0.414(3) 


0.626(7) 


0.01 


0.04 


0.04 


0.438(3) 


0.642(6) 


0.01 


0.05 


0.01 


0.387(3) 


0.613(9) 


0.01 


0.05 


0.02 


0.414(3) 


0.627(7) 


0.01 


0.05 


0.03 


0.440(3) 


0.642(6) 


0.01 


0.05 


0.04 


0.465(3) 


0.657(6) 


0.01 


0.05 


0.05 


0.489(3) 


0.672(5) 


0.02 


0.01 


0.01 


0.250(3) 


0.546(20) 


0.02 


0.02 


0.01 


0.292(3) 


0.560(14) 


0.02 


0.02 


0.02 


0.325(3) 


0.585(11) 


0.02 


0.03 


0.01 


0.327(3) 


0.579(12) 


0.02 


0.03 


0.02 


0.358(3) 


0.597(9) 


0.02 


0.03 


0.03 


0.386(3) 


0.615(8) 


0.02 


0.04 


0.01 


0.359(3) 


0.597(10) 


0.02 


0.04 


0.02 


0.388(3) 


0.618(8) 


0.02 


0.04 


0.03 


0.415(3) 


0.631(7) 


0.02 


0.04 


0.04 


0.440(3) 


0.649(6) 


0.02 


0.05 


0.01 


0.388(3) 


0.614(9) 


0.02 


0.05 


0.02 


0.415(3) 


0.631(7) 


0.02 


0.05 


0.03 


0.441(2) 


0.647(6) 


0.02 


0.05 


0.04 


0.466(2) 


0.663(6) 


0.02 


0.05 


0.05 


0.490(2) 


0.679(5) 


0.03 


0.01 


0.01 


0.251(3) 


0.599(24) 


0.03 


0.02 


0.01 


0.289(3) 


0.589(17) 


0.03 


0.02 


0.02 


0.325(3) 


0.613(12) 


0.03 


0.03 


0.01 


0.324(3) 


0.603(13) 


0.03 


0.03 


0.02 


0.355(3) 


0.616(10) 


0.03 


0.03 


0.03 


0.387(3) 


0.643(9) 


0.03 


0.04 


0.01 


0.356(3) 


0.618(11) 


0.03 


0.04 


0.02 


0.385(3) 


0.631(9) 


0.03 


0.04 


0.03 


0.412(3) 


0.655(8) 


0.03 


0.04 


0.04 


0.442(3) 


0.670(7) 


0.03 


0.05 


0.01 


0.386(3) 


0.633(10) 


0.03 


0.05 


0.02 


0.413(3) 


0.646(8) 


0.03 


0.05 


0.03 


0.438(3) 


0.660(7) 


0.03 


0.05 


0.04 


0.463(3) 


0.675(6) 


0.03 


0.05 


0.05 


0.488(2) 


0.690(6) 
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FIG. 6: Pseudoscalar effective masses from WL-WL, LL-LL 
and GL-GL correlation functions for the three ensembles. The 
valence quarks are degenerate, with the masses equal to the 
light sea quark mass of that ensemble. 



FIG. 7: Vector effective masses from WL-WL and GL-GL 
correlation functions for the three ensembles. The valence 
quarks are degenerate, with the masses equal to the light sea 
quark mass of that ensemble. 



tions of source/sink operators give masses which gener- 
ally agree, but there are a number of cases where the 
difference is somewhat outside of the error bars. Since 
we are concerned about effectively sampling phase space, 
the quoted statistical errors may be under-estimated due 
to possible long-range autocorrelations. We can use the 
fact that the different sources sample different parts of 
the lattice to compensate for this. For example, for the 
m sca _ q 02 ensemble, we find the mi = m 2 = 0.02 pseu- 
doscalar meson mass has a value of 0.3247(8), 0.3212(12) 
and 0.3251(18) from the PP channel, and 0.3254(12), 
0.3226(11) and 0.3287(19) from the AA channel. Taking 
these six values as independent and calculating a stan- 
dard deviation from them gives 0.3246(26), which has an 
error bar larger than from any of the individual mea- 
surements. From this example, amongst the ones with 
the largest differences, we have chosen to take the aver- 
ages of the results for the meson masses obtained from 
different measurements, wherever possible, as our best 
estimates of the meson masses. Since the different mea- 
surements involve different source locations on our lat- 
tices, the spread of values is also an indication of the 
statistical errors of our ensembles. We correspondingly 
inflate the largest individual error by a factor of 1.5 to 



bring it in rough agreement with this spread. The final 
results are summarized in Table |V1 where the values of 
771 p are taken from the averages of both the PP and AA 
channels in different measurements. These mass values 
will then be used in the subsequent discussions of the 
paper. 



Residual Mass 



The quantitative description of the explicit chiral sym- 
metry breaking for domain wall fermions is the residual 
@i EH 1 which is determined from the ratio 



mass, m rcs 



R(t) 



(Ea^OMMO))' 



(7) 



The definitions of Jg q (x,t) and Jg(x,t) for domain wall 
fermions can be found in Ref. [a]. We follow [ll| 
and determine m rcs at a given quark mass, denoted by 
TO res( m /)i by fitting R(t) to a constant over a time range 
where R(t) shows a good plateau. The results from the 
WL-WL correlation functions are shown in Table IVII 

While by definition m res should be a constant, the op- 
erator we choose to determine it happens to depend on 



10 
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FIG. 8: The residual mass determined from the WL-WL cor- 
relators. The solid line is a linear fit to the dynamical points 
with m va i = m; ca . 



the quark mass, as shown in Figure [8l This quark mass 
dependence comes from the higher-dimension operators 
in the divergence of the axial current, which should be 
suppressed by powers of a d, [111, HH • We choose to deal 
with these lattice artifacts by extrapolating m rcs (mj) to 
the zero quark mass limit in both the sea and valence 
sectors. Doing so gives 



mres = 0.00308(4). 



(8) 



This is about 1/3 of the lightest quark mass in our sim- 
ulations, which may seem larger than ideal. Neverthe- 
less, from a field theoretic point of view, to leading order 
the residual mass is just an additive renormalization to 
the input quark mass. Neglecting higher-order correc- 
tions, the chiral limit is defined at mj + Tires = or 
rrif = —m TCS , where is the input quark mass. This is 
the basis for the chiral extrapolations discussed in Sec- 
tion[V] In this paper we do not attempt to discuss higher- 
order corrections to the residual mass, which should be 
small in our simulations. 



D. Pseudoscalar Decay Constants 

The decay constant fp for a charged pseudoscalar me- 
son is defined by 



ifpq l >e-** = (0\A ll (x)\ir(q)), 



(9) 



where A^x) is the (partially) conserved axial vector cur- 
rent operator for domain wall fermions [8|, and 1 7r(qr) ) 



TABLE VI: Residual mass obtained with the Coulomb gauge 
fixed wall source. 



sea 






/„■ — f 

frnin ^max 


v 2 Ido f 


rn ' ( Tin e \ 


n ni 


01 


01 

U.Ul 


S 1 (S 
O 1U 




fl 0031 09(9^ 


n ni 
U.U1 


n no 
U.Uz 


n no 
U.Uz 


q i a 




U.UUzyDz^zo ) 


0.01 


0.03 


0.03 


8-16 


1.0(7) 




0.01 


0.04 


0.04 


8-16 


0.9(7) 


002752 flQl 


n no 


U.Ul 


U.Ul 




i r\(7\ 
L.yjy i ) 


U.UUOOOZ^ZO ) 




u.uz 


n n9 
u.uz 


C If! 


n q(7\ 


n no^i Q7f9n^ 


0.02 


0.03 


0.03 


8-16 


0.9(7) 


0.003070(19) 


0.02 


0.04 


0.04 


8-16 


0.9(7) 


0.002961(18) 


0.03 


0.01 


0.01 


8-16 


1.3(8) 


0.003479(27) 


0.03 


0.02 


0.02 


8-16 


1.5(9) 


0.003337(24) 


0.03 


0.03 


0.03 


8-16 


1.7(9) 


0.003201(22) 


0.03 


0.04 


0.04 


8-16 


1.9(10) 


0.003085(20) 



is the state vector for a pseudoscalar meson with mo- 
mentum q. We normalize the state vector such that 
Mg)K(g')) = {2n) 3 2E q S(q- q'). A^ is related to the 
local axial vector operator by A^ — ZaA^ for physics at 
long distances. It is conventional to determine the ma- 
trix element of (0|Ao|tt) from the amplitude of the LL-LL 
correlation functions [a [ll| , Ajf 1 Aq . Here a double-cosh 
fit to the LL-LL and GG-GG correlation functions was 
performed to get a more reliable estimate of the ampli- 
tude of the LL-LL correlation functions, fp is obtained 
by 



y mp 



(10) 



A novel way of determining fp is from the WL-WL 
correlators. In contrast to the point-like sources, the nor- 
malization of the Coulomb gauge fixed wall source is not 
known analytically, and we need to compute the relative 
normalization of the source operator to the conserved 
current numerically. This can be done by employing the 
following ratio, 



w 



(ii) 



where P w (0) is the Coulomb gauge fixed pseudoscalar 
source operator. The pseudoscalar decay constant can 
then be computed using the following formula: 



f (2) 

J p 



mp 



(12) 



This method has the advantage of giving smaller sta- 
tistical errors compared to the conventional way of de- 
termining fp from the point-like sources. We note that 
here N\y is merely a normalization factor for the two op- 
erators Ao(x,t) and Aq (x,t) computed under the same 
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TABLE VII: Z A from WL-WL correlation functions for all 
three ensembles. 





mi 


?7l2 


t i - tmax 


X 2 /dof 


Z A 


0.01 


0.01 


0.01 


8-16 


0.5(5) 


0.71807(14) 


0.01 


0.02 


0.02 


8-16 


7(6) 


0.71930(11) 


0.01 


0.03 


0.03 


8-16 


0.8(6) 


0.72071(9) 


0.01 


0.04 


0.04 


8-16 


0.9(7) 


0.72222(8) 


0.02 


0.01 


0.01 


8-16 


1.2(8) 


0.71799(14) 


0.02 


0.02 


0.02 


8-16 


1.0(7) 


0.71925(10) 


0.02 


0.03 


0.03 


8-16 


0.9(7) 


0.72069(8) 


0.02 


0.04 


0.04 


8-16 


0.8(7) 


0.72221(7) 


0.03 


0.01 


0.01 


8-16 


1.7(9) 


0.71806(15) 


0.03 


0.02 


0.02 


8-16 


2.0(10) 


0.71931(11) 


0.03 


0.03 


0.03 


8-16 


2.2(11) 


0.72073(10) 


0.03 


0.04 


0.04 


8-16 


2.2(10) 


0.72226(9) 


TABLE VIII: Pseudoscalar decay constants as 


computed from 


the double-cosh fits to the LL-LL and GG-GG axial vector 


correlation functions for all three ensembles. 




sea 

mi 




mi 


m.2 




Jp 


0.01 




0.01 


0.01 




0.0893(13) 


0.01 




0.01 


0.04 




0.1007(12) 


0.01 




0.04 


0.04 




0.1116(10) 


0.02 




0.02 


0.02 




0.1003(8) 


0.02 




0.02 


0.04 




0.1070(6) 


0.02 




0.04 


0.04 




0.1133(7) 


0.03 




0.03 


0.03 




0.1094(13) 


0.03 




0.03 


0.04 




0.1111(12) 


0.04 




0.04 


0.04 




0.1143(12) 



condition. Unlike Za, which is a constant up to correc- 
tions of 0(a 2 ), the mass dependence of Nw has physical 
significance and thus it is the value of N\y evaluated at 
each valence quark mass that should be used in Eq . (fT2")) . 

The value of Za is determined using the method de- 
scribed in [ll[ and [1]. We quote the results obtained 
from the WL-WL functions in Table IVIII Results from 
other source-sink combinations are in good agreement 
with the shown values, although we do not display them 
in the paper. Taking the chiral limit removes an 0(a 2 ) 
lattice artifact and gives 

Z A = 0.7162(2). (13) 

This is then used to extract fp from the axial vector 
correlators using either Eq. (fT0]) or Eq. (fl2"|) , the results 
of which are given in Tables I Villi and IIXI respectively. 
These two methods of extracting fp agree well, with the 
errors on fp being smaller than fp^ in most cases. 



TABLE IX: Pseudoscalar decay constants as computed from 
the WL-WL axial vector correlation functions for all three 
ensembles. 





mi 


7712 


r(2) 
■IP 


0.01 


0.01 


0.01 


0.0895(7) 


0.01 


0.02 


0.02 


0.0977(7) 


0.01 


0.03 


0.03 


0.1046(7) 


0.01 


0.04 


0.04 


0.1108(7) 


0.02 


0.01 


0.01 


0.0939(7) 


0.02 


0.02 


0.02 


0.1006(6) 


0.02 


0.03 


0.03 


0.1067(6) 


0.02 


0.04 


0.04 


0.1123(6) 


0.03 


0.01 


0.01 


0.0980(6) 


0.03 


0.02 


0.02 


0.1042(6) 


0.03 


0.03 


0.03 


0.1101(6) 


0.03 


0.04 


0.04 


0.1157(6) 



V. OBSERVABLES IN THE CHIRAL LIMIT 

Since we perform our simulations at quark masses that 
are heavier than the physical up and down quark masses, 
extrapolations are needed to obtain physical values for 
various quantities. Chiral perturbation theory is the cor- 
rect effective theory to describe low-energy QCD and 
should be used to guide the extrapolations. However, the 
quark masses in our simulations are too heavy for next- 
to-leading-order (NLO) chiral perturbation theory to be 
consistent with our results |25j. Although independent 
NLO fits for the pseudoscalar meson masses and decay 
constants are consistent with our data, our attempt to fit 
both simultaneously fails badly. This is to be expected 
since the pion masses at the light dynamical points in our 
simulations range from 400 MeV to 627 MeV. This mass 
scale is likely to be outside of the region where chiral per- 
turbation theory is applicable. As will be demonstrated 
in the following, our data exhibits linear behavior in the 
range of quark masses in our simulations. Thus we resort 
to simple linear extrapolations to obtain physical results 
in the light quark mass limit. In many cases this co- 
incides with the predictions of the leading-order chiral 
perturbation theory, but in fact the reason we are doing 
this is based purely on phenomenological observations, 
not from the underlying effective theory. 



A. Determination of the Lattice Scale 

We start out by setting the lattice scale for these simu- 
lations from the vector p mass. Although at the lightest 
quark masses the vector meson mass in our simulations is 
slightly larger than the 2m w threshold, the requirement 
for non-zero relative momentum for the two pions pro- 
hibits the vector meson from decaying. We have chosen 
to do a partially quenched linear fit to the values of my in 
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FIG. 9: Linear extrapolations of the vector mass using 
Eq. (114[) . Here the circles are partially quenched data points 
from the m; ca = 0.01 ensemble, the squares are the unitary 
points with mi = mi = m; ca and the diamond is the value of 
mv in the chiral limit. 



Table fVl restricted to mi (7712) < 0.03, with the following 
phcnomcnological form : 

m v = A {mf a + m TCS ) + B (mi + m 2 + 2m ros ) + C. (14) 

The values of the fit parameters are shown in Table [X] 
As an illustration, we show the fit curves through the 

m? e 



= 0.01 data points [29[ and the unitary points in 
Figure [H Setting mi = m 2 = m^ ca = — m res gives the 
p mass in the chiral limit, from which we determine the 
lattice scale to be 



1.61(3) GeV. 



(15) 



We also determined the lattice scale from the static 
quark potential using the Coulomb gauge method [26| . 
The preliminary analysis of these ensembles was reported 
in [27| . More detailed analysis of this calculation is the 
subject of another paper in progress and will not be ad- 
dressed here. We simply quote the updated value of the 
lattice scale determined using rg = 0.5 fm, which is 



1.63(5) GeV. 



(16) 



Our third choice for determining the lattice scale is 
the "method of lattice planes" [28| . In this approach the 
vector meson mass is plotted against the pseudoscalar 
meson mass squared and the interception of this data 
with the physical point (m^ , mjf,) is found. Since the 
lattice data straddles the kaon data, a chiral interpolation 



FIG. 10: Linear extrapolation of a 1 versus m; ca from the 
method of planes. 



rather than extrapolation is required in the valence quark 
sector. The results for the inverse lattice spacing for each 
of the sea quark masses is shown in Figure 1101 Shown 
in this figure is a linear extrapolation in the sea quark 
mass to the chiral limit, mp a = — m rcs . In this limit, we 
obtain 



Q (lattice plane 



1.62(5) GeV. 



(17) 



All these methods give consistent results for the lattice 
scale. However, each method has its own limitations. 
The physical p meson is unstable and has the probability 
of decaying into two pions (though it is not the case in 
our simulations). The static quark potential suffers a 
small (few percent) uncertainty concerning the value of 
r . For the method of lattice planes, the extrapolation in 
the sea sector is purely phenomenological. Nevertheless, 
the three determinations of a" 1 agree well, so we take 
the average of these for our central value, and take as an 
error an average of the statistical errors. This gives 



1.62(4) GeV, 



(18) 



which will be used whenever a lattice scale is needed, and 
the errors will be propagated accordingly by quadrature. 



B. Quark Masses and Decay Constants 

A precise determination of the physical quark masses 
demands well-controlled chiral extrapolations in the light 
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FIG. 11: Linear extrapolation of the pseudoscalar meson mass 
through the light dynamical points. The crosses are the phys- 
ical 7r and K points from which we obtain the values of rh and 
m 3 . 



FIG. 12: Linear extrapolation of the pseudoscalar meson de- 
cay constant through the light dynamical points. The dia- 

(2) 

monds are the values of f n and fx obtained from p p ' . 



quark limit. As already mentioned, the relatively heavy 
quark masses in our simulations prevent us from using 
the next-to-leading order chiral perturbation theory to do 
the extrapolations. Again we restrict ourselves to linear 
extrapolation. We chose to fit only to the dynamical 
points where m va \ = rof ea , using the following formula: 

m% = 2B(mf a + m res ). (19) 

The value of B is given in Table |XJ We then used the 
physical pion and kaon masses as inputs to determine the 
average light quark mass and the strange quark mass, as 
shown in Figure 1111 The bare quark masses in lattice 
units are listed in the following: 

rh = 0.00162(8) (20) 
m s = 0.0390(21). (21) 

The above results have combined the bare input quark 
masses and the residual mass, and the quoted errors come 
from the errors on the pseudoscalar meson masses, the 
residual mass, and the lattice scale. Systematic uncer- 
tainties associated with the finite lattice spacing, chiral 
extrapolation and e-m splitting, etc. have not been inves- 
tigated yet. We will employ the non-perturbative renor- 
malization (NPR) technique Q to obtain the renormal- 
ized quark masses. The NPR analysis is still in progress. 

To obtain the physical results of f„ and fx, we fit the 
results of the pseudoscalar decay constants in Tables IIXI 
and I Villi to the following linear form separately: 

/ P = /o + L(mP + m te8 ). (22) 



We can see from Figure [12] that the data points are quite 
linear with our current statistics. The fitting parameters 
are also tabulated in Table [X] Using the quark masses 
and lattice scale determined in the above, we found that 
the physical values of f n , fx and fx/ to be 126.9(45) 
MeV, 157.3(51) MeV and 1.240(31) respectively using 

fp\ Similar analysis using /J?^ gives /„■ = 126.7(36) 
MeV, fx = 157.6(45) MeV and f K /U = 1.244(20), 
which are consistent with the results from fp . We take 
the averages of these two methods as our best estimates 
for these quantities, which are summarized in the follow- 
ing: 

U = 127(4) MeV, (23) 
fx = 157(5) MeV, (24) 
JkIU = 1.24(2). (25) 

Our convention for the pseudoscalar decay constant is 
such that the experimental value of /„. is about 131 MeV, 
f K about 160 MeV and f K /U about 1.22. Thus our 
results are in good agreement with the experiment. 
We stress that these physical results are obtained in 



TABLE X: Linear fit parameters for various quantities. 



my 


A 


= 1.42(23) 


B = 


= 1.58(16) 


C = 0.483(9) 


mp 






B = 


2.285(24) 




f ( p 1} 


L 


= 1.01(9) 


fo = 


0.0767(20) 




f(2) 
JP 


L 


= 1.02(5) 


fo = 


0.0765(11) 





14 



the context of simple linear extrapolations. Systematic 
uncertainties may arise from the presence of chiral logs. 
More sophisticated examination of the chiral behavior of 
various observables requires simulations at lighter quark 
masses, which will be further addressed in the near fu- 
ture. Nevertheless, such agreement is encouraging, and 
demonstrates that simulations with 2+1 flavors of do- 
main wall fermions are promising. 

VI. COMPARISON OF RHMC I AND RHMC II 

As discussed in Section ing there seem to be some long- 
range autocorrelations in these ensembles generated with 
2+1 flavors of domain wall fermions. However, early 
simulations with two flavors of dynamical domain wall 
fermions did not have this problem 22]. The biggest 
difference between those two-flavor simulations and this 
work lies in the way the following ratio is evaluated, 



det 


D dwf ( M 5 , m s ) L>dwf (M 5 , m s ) 


det 


^dwf( M 5,1)#dwf(M 5 ,1) 





In the two-flavor case, this ratio was evaluated using one 
pseudo-fermion field in the HMC evolution [22j, which 
was found to reduce the stochastic noise substantially, 
and a larger step size could be used. At the beginning 
of these simulations this method had not yet been im- 
plemented in the RHMC algorithm, and two separate 
pseudo-fermion fields were needed to evaluate Eq. (|26|) . 
Upon completion of the work reported in the previous 
sections, the new code was ready, and we were motivated 
to perform an extension of the m^ 63, = 0.03 ensemble us- 
ing this new code to investigate any effects it might have 
on the ensemble properties. This variant of the RHMC 
algorithm, RHMC II, has been described in detail in Sec- 
tion HI From Table Q] we can see that the coarsest step 
size, St, of RHMC II is almost twice as large as RHMC 
I, and fewer steps of intermediate and finest integration 
sizes are needed. We even obtained higher acceptance 
with these parameters. There is clearly at least a factor 
of two speedup by switching from RHMC I to RHMC II. 

Figure [T"3l shows the comparison of the evolution of the 
topological charge on the lattices generated using RHMC 
I and RHMC II. The tunneling of topological charge in 
RHMC II appears to be much better than RHMC I. The 
long-range autocorrelations present in RHMC I do not 
show up in RHMC II. This can be partly attributed to 
the use of a larger step size in the molecular dynamics 
evolution. Intuitively, a larger step size allows more effi- 
cient movement through, and the evolution is less likely 
to be trapped in a small region of, phase space. This con- 
sequently helps to move the topological charge from one 
sector to another. However, the histograms of the topo- 
logical charge from these two segments of the ensemble 
are quite similar, indicating that RHMC I, while show- 
ing longer range autocorrelations, does give the same re- 
sults as the improved version of the algorithm, RHMC 




1000 2000 3000 4000 5000 6000 7000 8000 
MD time unit 




top ^top 



FIG. 13: The comparison of the topological charge evolution 
and corresponding histogram measured from the two portions 
of the m; ea = 0.03 ensemble, in which different algorithms 
(RHMC I and RHMC II) were used to evaluate Eq.JT]). The 
dashed line indicates where the algorithm was changed. 



II. We have also seen large fluctuations of the effective 
masses in Section |TV] As a comparison, we show the ef- 
fective masses for the pseudoscalar and vector WL-WL 
correlation functions, measured every 20 trajectories with 
t src = at the unitary point in Figure [14] While the vec- 
tor effective mass remains noisy for RHMC II, there is no- 
ticeable improvement in the pseudoscalar effective mass, 
where a better plateau can be identified from t = 8 to 
t = 24. 

We believe that RHMC II not only speeds up the gauge 
field generation by more than a factor of 2, but also helps 
to move through phase space more quickly and improve 
the quality of the ensembles we generate. The interplay 
of these two factors makes it a more favorable algorithm 
to use than RHMC I. 



VII. CONCLUSIONS 

We have presented the results of light meson physics 
from lattice simulations with 2+1 flavors of dynamical 
domain wall fermions at a lattice spacing of about 0.12 
fm in a (2 fm) 3 volume. We used the improved RHMC 
algorithm to generate the gauge configurations, which 
turned out to be quite efficient in reducing the comput- 
ing cost. Although there are some long-range autocorre- 
lations in these ensembles, we have been able to obtain 
physical results with better than 5% statistical accuracy. 
The residual chiral symmetry breaking in these calcula- 
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FIG. 14: The comparison of the effective masses for the pseu- 
doscalar and vector WL-WL correlation functions at the uni- 
tary point from the two portions of the m; ea = 0.03 ensemble 
(RHMC I and RHMC II). 



could result in sizable finite volume effects. To pursue 
this issue further, the RBC and UKQCD collaborations 
are performing 2+1 flavor simulations on a 24 3 x 64 vol- 
ume, or about (3 fm) 3 in physical units, at about the 
same lattice spacing. This larger volume allows calcu- 
lations with lighter quark masses, from which we hope 
to investigate the light quark limit in the context of chi- 
ral perturbation theory. Comparing the physical results 
in two different volumes will also enable us to quantita- 
tively study finite volume effects. This will be especially 
important for baryon physics. Secondly, more than one 
lattice spacing is needed to investigate the scaling be- 
havior for domain wall fermions and the Iwasaki gauge 
action, and the continuum limit should be taken to re- 
move lattice artifacts. We will address this question in 
upcoming simulations in a (3 fm) 3 volume at a weaker 
coupling. 

With good control over chiral symmetry and an ex- 
act and efficient algorithm, we are vigorously pursuing 
2+1 flavor domain wall QCD. The theoretical promise of 
this approach is now producing numerical results for the 
physically interesting 2+1 flavor case. 



Acknowledgments 



tions is small compared to the earlier calculation with 
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APPENDIX A: TABLES OF MESON MASSES 



TABLE XL Pseudoscalar meson masses from GL-GL correla- 
tion functions as determined from the pseudoscalar channel. 





mi 


m 2 tm 


n ^max 


X 2 /dof 


mp 


0.01 


0.01 


0.01 


8-16 


0.6(6) 


0.2460(10) 


0.01 


0.01 


0.04 


8-16 


0.8(6) 


0.3545(8) 


0.01 


0.04 


0.04 


8-16 


1.1(7) 


0.4382(7) 


0.02 


0.02 


0.02 


8-16 


0.9(7) 


0.3247(8) 


0.02 


0.02 


0.04 


8-16 


0.9(7) 


0.3868(8) 


0.02 


0.04 


0.04 


8-16 


0.9(7) 


0.4410(7) 


0.03 


0.03 


0.03 


8-16 


3.8(15) 


0.3907(15) 


0.03 


0.03 


0.04 


8-16 


3.6(14) 


0.4184(14) 


0.03 


0.04 


0.04 


8-16 


3.5(14) 


0.4446(13) 


TABLE XII: 


Pseudoscalar 


meson masses from GL-GL corre- 


lation functions as determined from the axial vector channel. 


m s r 


mi 


Tfl2 ^min ^max 


x 2 /dof 


nip 


0.01 


0.01 


0.01 


8-16 


2.4(10) 


0.2467(15) 


0.01 


0.01 


0.04 


8-16 


1.9(8) 


0.3547(14) 


0.01 


0.04 


0.04 


8-16 


2.0(9) 


0.4379(10) 


0.02 


0.02 


0.02 


9-16 


2.1(10) 


0.3254(12) 


0.02 


0.02 


0.04 


9-16 


1.8(9) 


0.3877(10) 


0.02 


0.04 


0.04 


9-16 


1.7(9) 


0.4413(10) 


0.03 


0.03 


0.03 


7-16 


3.2(11) 


0.3881(16) 


0.03 


0.03 


0.04 


7-16 


2.8(10) 


0.4170(16) 


0.03 


0.04 


0.04 


7-16 


1.7(9) 


0.4425(13) 



TABLE XIII: Vector meson masses from GL-GL correlation 
functions for all three ensembles. 





mi 


m-2 


tmin 'max 


X 2 /dof 


my 


0.01 


0.01 


0.01 


10-16 


0.7(7) 


0.558(12) 


0.01 


0.01 


0.04 


10-16 


0.6(6) 


0.6016(45) 


0.01 


0.04 


0.04 


10-16 


0.7(7) 


0.6443(23) 


0.02 


0.02 


0.02 


10-16 


1.4(10) 


0.5914(69) 


0.02 


0.02 


0.04 


10-16 


0.7(8) 


0.6217(47) 


0.02 


0.04 


0.04 


10-16 


0.4(6) 


0.6524(35) 


0.03 


0.03 


0.03 


10-16 


0.7(8) 


0.6519(75) 


0.03 


0.03 


0.04 


10-16 


0.7(7) 


0.6643(63) 


0.03 


0.04 


0.04 


10-16 


0.6(7) 


0.6771(54) 
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TABLE XIV: Pseudoscalar meson masses as determined from TABLE XVII: Pseudoscalar meson masses from WL-WL cor- 
the WL-WL pseudoscalar correlation functions. relation functions with P+A boundary conditions as deter- 



mf a 


TTii 
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k m in v ill ax 
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mp 


mined from the pseudoscalar channel. 






0.01 


0.01 
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1 1(9) 


0.2474(16) 


mf a 


rrii 


nio 
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n m 


o^i 

U.UO 


09 

U.UZ 


1 99 

1U~ZZ 


0.4161 


18) 














U.Ul 


O^i 
u.uu 


0*3 

U.UO 


1 99 


0.4418 


18) 














01 

U.Ul 


O^i 

U.UO 


04 

U.U4: 


1 99 

1U~ZZ 


0.4666 


17) 


TABLE XV: Pseudoscalar meson masses from WL-WL corre- 


01 

U.Ul 


0^ 

U.UO 


0^ 

U.UO 


10-22 


0.4906 


17) 


lation functions as determined from the axial vector channel. 


09 
u.uz 


01 
U.Ul 


01 
U.Ul 


1 99 
iu-zz 


0.2527 


20) 


m B r 


mi 


mi 


^min Cmax 


x 2 /dof 


rap 


09 
u.uz 


09 

U.UZ 


01 

U.Ul 


1 99 

1U~ZZ 


0.2910 


19) 


0.01 


0.01 


0.01 


8-16 


0.8(7) 


0.2466(16) 


09 

u.uz 


09 
U.UZ 


09 
U.UZ 


1 99 
iu-zz 


0.3251 


18) 


0.01 


0.02 


0.02 


8-16 


0.6(6) 


0.3215(14) 


09 
U.UZ 


0*3 
U.Uo 


01 
U.Ul 


1 99 
iu-zz 


0.3253 


19) 


0.01 


0.03 


0.03 


8-16 


0.8(7) 


0.3828(13) 


09 

U.UZ 


0*3 

U.UO 


09 

U.UZ 


10-22 


0.3564 


17) 


0.01 


0.04 


0.04 


8-16 


1.2(8) 


0.4365(12) 


0.02 


0.03 


0.03 


10-22 


0.3855 


17) 


0.02 


0.01 


0.01 


8-16 


1.4(9) 


0.2479(13) 


09 
u.uz 


04 

U.U4: 


01 

U.Ul 


1 99 

1U~ZZ 


0.3567 


19) 


0.02 


0.02 


0.02 


8-16 


1.6(10) 


0.3226(11) 


O 09 
U.UZ 


n 0/i 

U.U4 


n 09 
u.uz 


i n 99 

lU-zz 


0.3856 


17) 


0.02 


0.03 


0.03 


8-16 


1.7(10) 


0.3842(10) 


09 
u.uz 


04 

U.U4: 


0*3 

U.UO 


1 99 

1U~ZZ 


0.4129 


16) 


0.02 


0.04 


0.04 


8-16 


1.9(10) 


0.4385(9) 


0.02 


0.04 


0.04 


10-22 


0.4390 


151 


0.03 


0.01 


0.01 


8-16 


2.2(11) 


0.2526(14) 


0.02 


0.05 


0.01 


10-22 


0.3859 


19) 


0.03 


0.02 


0.02 


8-16 


3.0(13) 


0.3268(14) 


n no 
U.Uz 


ft (If! 

U.Uo 


n no 
U.Uz 


i n oo 
lU-zz 


0.4131 


17) 


0.03 


0.03 


0.03 


8-16 


2.9(13) 


0.3884(13) 


0.02 


0.05 


0.03 


10-22 


0.4390 


16) 


0.03 


0.04 


0.04 


8-16 


2.8(13) 


0.4427(11) 


0.02 


0.05 


0.04 


10-22 


0.4639 


15) 














0.02 


0.05 


0.05 


10-22 


0.4878 


15) 














0.03 


0.01 


0.01 


10-22 


0.2496 


21) 














0.03 


0.02 


0.01 


10-22 


0.2885 


19) 


TABLE XVI: 


Vector meson masses from WL-WL correlation 


0.03 


0.02 


0.02 


10-22 


0.3227 


19) 


functions for all three ensembles. 






0.03 


0.03 


0.01 


10-22 


0.3235 










X 2 /dof 




19) 


m s r 


mi 


mi 


^min ^max 


my 


0.03 


0.03 


0.02 


10-22 


0.3544 


18) 


0.01 


0.01 


0.01 


9-16 


1.2(9) 


0.542(11) 


0.03 


0.03 


0.03 


10-22 


0.3837 


18) 


0.01 


0.02 


0.02 


9-16 


1.9(11) 


0.5742(58) 


0.03 


0.04 


0.01 


10-22 


0.3556 


19) 


0.01 


0.03 


0.03 


9-16 


1.9(11) 


0.6078(38) 


0.03 


0.04 


0.02 


10-22 


0.3841 


18) 


0.01 


0.04 


0.04 


9-16 


1.6(10) 


0.6402(29) 


0.03 


0.04 


0.03 


10-22 


0.4116 


17) 


0.02 


0.01 


0.01 


6-16 


0.8(6) 


0.5533(44) 


0.03 


0.04 


0.04 


10-22 


0.4379 


17) 


0.02 


0.02 


0.02 


6-16 


0.8(6) 


0.5845(30) 


0.03 


0.05 


0.01 


10-22 


0.3854 


19) 


0.02 


0.03 


0.03 


6-16 


0.9(6) 


0.6162(23) 


0.03 


0.05 


0.02 


10-22 


0.4120 


18) 


0.02 


0.04 


0.04 


6-16 


1.1(7) 


0.6475(19) 


0.03 


0.05 


0.03 


10-22 


0.4380 


17) 


0.03 


0.01 


0.01 


10-16 


0.9(8) 


0.621(16) 


0.03 


0.05 


0.04 


10-22 


0.4631 


17) 


0.03 


0.02 


0.02 


10-16 


0.8(8) 


0.6253(82) 


0.03 


0.05 


0.05 


10-22 


0.4873 


16) 


0.03 


0.03 


0.03 


10-16 


1.0(9) 


0.6464(55) 














0.03 


0.04 


0.04 


10-16 


1.4(11) 


0.6718(42) 
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TABLE XVIII: Pscudoscalar meson masses from WL-WL cor- 
relation functions with P+A boundary conditions as deter- 
mined from the axial vector channel. 





mi 


m,2 


train ~ tmax 


mp 


0.01 


0.01 


0.01 


10-22 


0.2473(23) 


0.01 


0.02 


0.01 


10-22 


0.2878(21) 


0.01 


0.02 


0.02 


10-22 


c\ 0001 /i a\ 

0.3231(19) 


0.01 


0.03 


0.01 


10-22 


0.3236(20) 


0.01 


0.03 


0.02 


10-22 


0.3552(18) 


0.01 


0.03 


0.03 


10-22 


0.3847(17) 


0.01 


0.04 


0.01 


10-22 


0.3560(20) 


0.01 


0.04 


0.02 


10-22 


0.3850(17) 


n m 
0.01 


n n a 
U.04 


n no 
U.U3 


i n oo 
W-ZZ 


a a i ( ~\ c\ 
U.4124(loj 


0.01 


0.04 


0.04 


10-22 


0.4386(15) 


n m 


ft (If; 

U.Uo 


n m 
U.U1 


i n oo 
W-ZZ 


U.oo6U(19 ) 


n m 
0.01 


n n c 
U.Uo 


n no 
U.U2 


i n oo 
W-ZZ 


n /i i on/ 1 ^ 
U.41oU(l / ) 


n ni 
0.01 


n n c 
U.Oo 


n no 
U.U3 


i n oo 
W-ZZ 


U. 4,588(16) 


n m 
U.U1 


ft ftf! 

U.Uo 


n n a 
U.U4 


i n oo 


U.40O / (lo) 


n m 
0.01 


n n c 
U.Uo 


n n k 
U.Uo 


i n oo 
1U-22 


n /I Q77/1 c 

U.48/ /(lo) 


0.02 


0.01 


0.01 


10-22 


0.2535(22) 


0.02 


0.02 


0.01 


10-22 


0.2938(20) 


0.02 


0.02 


0.02 


10-22 


0.3287(19) 


0.02 


0.03 


0.01 


10-22 


0.3292(20) 


0.02 


0.03 


0.02 


10-22 


0.3605(18) 


0.02 


0.03 


0.03 


10-22 


0.3897(17) 


0.02 


0.04 


0.01 


10-22 


a oai a/i a\ 

0.3612(19) 


0.02 


0.04 


0.02 


10-22 


0.3899(17) 


n no 
0.02 


n n a 
U.U4 


n no 
U.U3 


i n oo 
W-ZZ 


A A 1 71 /l a\ 

U.41 / 1(10) 


n no 
U.Uz 


n n a 
U.U4 


n n a 
U.U4 


i n oo 


r\ a a on/i t:\ 
U. 4401^10 ) 


0.02 


0.05 


0.01 


10-22 


a oaa7/'ia'\ 

0.3907(19) 


n no 
0.02 


n n k 
U.Uo 


n no 
U.U2 


i n oo 
1U-22 


U.41 /o(l / ) 


n no 
0.02 


n n c 
U.Uo 


n no 
U.U3 


i n oo 
1U-22 


n /i/ioi/ic;\ 
U. 44,51(16) 


n no 
0.02 


n n k 
U.Uo 


n n /i 
U.U4 


i n oo 
1U-22 


U.46/ /(lo) 


n no 
0.02 


n n c 
U.Uo 


n n c 
U.Uo 


i n oo 
1U-22 


U. 4914(15) 


0.03 


0.01 


0.01 


10-22 


0.2486(23) 


0.03 


0.02 


0.01 


10-22 


0.2886(22) 


0.03 


0.02 


0.02 


10-22 


0.3235(21) 


0.03 


0.03 


0.01 


10-22 


0.3242(22) 


0.03 


0.03 


0.02 


10-22 


0.3555(20) 


0.03 


0.03 


0.03 


10-22 


0.3849(19) 


0.03 


0.04 


0.01 


10-22 


0.3565(22) 


0.03 


0.04 


0.02 


10-22 


0.3852(20) 


0.03 


0.04 


0.03 


10-22 


0.4126(19) 


0.03 


0.04 


0.04 


10-22 


0.4388(18) 


0.03 


0.05 


0.01 


10-22 


0.3863(22) 


0.03 


0.05 


0.02 


10-22 


0.4131(20) 


0.03 


0.05 


0.03 


10-22 


0.4389(19) 


0.03 


0.05 


0.04 


10-22 


0.4637(18) 


0.03 


0.05 


0.05 


10-22 


0.4877(17) 



TABLE XIX: Vector meson masses from WL-WL correlation 
functions with P+A boundary conditions for all three ensem- 



bles. 


m; ea 


mi 


m2 


£min - tmax 


mv 


0.01 


0.01 


0.01 


6-14 


0.546(13) 


0.01 


0.02 


0.01 


6-14 


0.563(9) 


0.01 


0.02 


0.02 


6-14 


0.579(7) 


0.01 


0.03 


0.01 


6-14 


0.580(7) 


0.01 


0.03 


0.02 


6-14 


0.595(6) 


0.01 


0.03 


0.03 


6-14 


0.610(5) 


0.01 


0.04 


0.01 


6-14 


0.597(6) 


0.01 


0.04 


0.02 


6-14 


0.611(5) 


0.01 


0.04 


0.03 


6-14 


0.626(5) 


0.01 


0.04 


0.04 


6-14 


0.641(4) 


0.01 


0.05 


0.01 


6-14 


0.613(6) 


0.01 


0.05 


0.02 


6-14 


0.627(5) 


0.01 


0.05 


0.03 


6-14 


0.642(4) 


0.01 


0.05 


0.04 


6-14 


0.657(4) 


0.01 


0.05 


0.05 


6-14 


0.672(3) 


0.02 


0.01 


0.01 


6-14 


0.538(13) 


0.02 


0.02 


0.01 


6-14 


0.560(10) 


0.02 


0.02 


0.02 


6-14 


0.579(7) 


0.02 


0.03 


0.01 


6-14 


0.579(8) 


0.02 


0.03 


0.02 


6-14 


0.597(6) 


0.02 


0.03 


0.03 


6-14 


0.614(5) 


0.02 


0.04 


0.01 


6-14 


0.597(7) 


0.02 


0.04 


0.02 


6-14 


0.614(5) 


0.02 


0.04 


0.03 


6-14 


0.631(5) 


0.02 


0.04 


0.04 


6-14 


0.647(4) 


0.02 


0.05 


0.01 


6-14 


0.614(6) 


0.02 


0.05 


0.02 


6-14 


0.631(5) 


0.02 


0.05 


0.03 


6-14 


0.647(4) 


0.02 


0.05 


0.04 


6-14 


0.663(4) 


0.02 


0.05 


0.05 


6-14 


0.679(4) 


0.03 


0.01 


0.01 


6-14 


0.577(15) 


0.03 


0.02 


0.01 


6-14 


0.589(11) 


0.03 


0.02 


0.02 


6-14 


0.601(8) 


0.03 


0.03 


0.01 


6-14 


0.603(9) 


0.03 


0.03 


0.02 


6-14 


0.616(7) 


0.03 


0.03 


0.03 


6-14 


0.630(6) 


0.03 


0.04 


0.01 


6-14 


0.618(7) 


0.03 


0.04 


0.02 


6-14 


0.631(6) 


0.03 


0.04 


0.03 


6-14 


0.645(5) 


0.03 


0.04 


0.04 


6-14 


0.660(4) 


0.03 


0.05 


0.01 


6-14 


0.633(6) 


0.03 


0.05 


0.02 


6-14 


0.646(5) 


0.03 


0.05 


0.03 


6-14 


0.660(5) 


0.03 


0.05 


0.04 


6-14 


0.675(4) 


0.03 


0.05 


0.05 


6-14 


0.690(4) 



